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The dynamics of critical slope self-organized critical models is studied, using a previous mapping 
into a linear interface depinning model dragged at one end. The model is solved obtaining the 
complete set of scaling exponents. Some results are supported by previous RG developed for constant 
force linear interface depinning models but others, like the linear dependency of the susceptibility 
with system size, are intrinsic of this model which belongs to a different universality class. The 
comparison of our results with numerical simulations of ricepile and vortexpile models reported in the 
literature reveals that, as in the constant force case, there are two universality classes corresponding 
to random and periodic pinning. 
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I. INTRODUCTION 



There are many natural phenomena where the dynam- 
ics does not take place smoothly but in the form of bursts 
or avalanches. Avalanche dynamics has been observed in 
earthquakes dynamics Q , granular piles , supercon- 
ductor vortex piles Q, the Barkhausen effect [||, crack 
propagation |7j, and more. The avalanches are char- 
acterized by their size s and duration T which in gen- 
eral follow the power law distributions P(s) ~ s~ Ts and 
p(T) ~ T- T \ respectively. 

Bak, Tang, and Wiesenfeld || introduced the notion 
of self-organized criticality (SOC) to explain this kind of 
dynamics, which seems to be common to many "differ- 
ent" phenomena. In their early formulation the critical 
state of a SOC system is an attractor of its dynamics 
and therefore there is no need of fine-tuning to reach the 
critical state. This idea was illustrated with very simple 
cellular automata, the sandpile models. In sandpile mod- 
els an integer p|-|l^] or continuous 13 - TBI variable Zi is 
defined in a d-dimensional lattice. Zi can be the height 
of grains columns [8|]9|, |i^ , |l6| , the stress accumulated in 
certain fault JTJJ], the vortex density etc. A critical 
threshold condition for local relaxation (toppling) is also 
given, which can be a critical height Zi > z c P- 5^ , H , 
critical slope m; = 2, — 2^+1 > m c |BO,E5l, or critical 

Laplacian U = J2nn Zi ~ - > ^ c ^ nresn0 ^ con ~ 
dition. When this condition is fulfilled the site relaxes 
transferring grains to its nearest neighbors (nn), other- 
wise certain amount of grains is added from an external 
field. 

Critical height models have been extensively studied in 
the literature, either by mean field ||18|1 a nd field theories 
fL9f , renormalization group (RG) |20| , |21| , and numerical 
simulations |p|- ]l0| , ^3| , ^4[ . However, the study of critical 
slope models has been in general limited to numerical 
simulations [^11 lq ], which are not so extensive as in 
the critical height variant, while theoretical approaches 
are almost absent. There are many experimental situa- 
tions, sandpiles and vortexpiles for instance, where the 



critical slope threshold condition is more appropriate, 
requiring a rigorous treatment which captures the gen- 
eral features of critical slope models. In this direction 
Paczuski and Boettcher |2^] provided an starting point. 
They mapped a critical slope model into a linear interface 
depinning (LID) model where the interface is dragged at 
one end, showing the existence of common behavior be- 
tween these models. The importance of their observation 
is that LID models can be solved using continuous ap- 
proaches |23|-^5||. Their work was inconclusive because, 
up to our knowledge, the problem of an interface moving 
on a random environment and dragged at one end has not 
been solved. However, they observed that some scaling 
exponents, but not all, are identical to those measured 
for an interface driven uniformly. 

In the present work we solve the problem of an inter- 
face moving on a random environment and dragged at 
one end. We start with a LID equation with the appro- 
priate boundary conditions. We then look for the fluctu- 
ations around the average in the stationary state. In the 
thermodynamic limit these fluctuations satisfy an equa- 
tion similar to that of an interface driven uniformly. The 
critical exponents are then computed using previous RG 
calculations for the uniformly driven interface [^3] p5| . 
The dynamic and roughness scaling exponents are found 
identical while other exponents become different. These 
results show a very good agreement when compared with 
numerical estimates of the scaling exponents reported in 
the literature. 



II. THE MODEL 

Consider the following one-dimensional critical slope 
model. An integer variable Zi (i = 1, 2, . . . , L) is defined. 
A site where the slope rrii — Zi — z i+ i is greater than a 
threshold m c is said to be active and relaxes according 
to the toppling rule 



Z{ ► z>i 1 , 
Zi+i -> z i+ i + 1. 



(1) 
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Grains are added at site i = at rate c while the other 
boundary i = L — 1 is open. Usually c is assumed very 
small in such a way that a new grain is added only when 
no site is active, which is the usual separation of time 
scales assumed in numerical simulations. 

Paczuski and Boettcher j2^] noted that instead of fol- 
low the evolution of the column heights zi one can follow 
the evolution of hi(t), the total number of toppling events 
at site i up to time t. The evolution rules for hi can be 
easily obtained from those for Zi . If we start with an ini- 
tial configuration where Zj = the number of grains at 
site i is obtained adding the number of grains received 
from the left neighbor (hi-\) and subtracting the number 
of grains transferred to the right neighbor (hi), i.e. 

Zi = h,-i - hi. (2) 

Using this expression one may compute the local slope 
mi — Zi — Zi + i, resulting 

m< - m c — hi+i + hi-i - 2h, - m c , (3) 

On the other hand, hi increases by one unit every time 
the site i topples, i.e. every time m t > m c , which can be 
written as 

hi(t + l)-hi(t) = Q(mi-m e ), (4) 

where Q(x) is the Heaviside unit step function. 

The evolution rules in eqs. ([|) and dJ) are just the 
discretized evolution rules for LID models |2(|, taking 
m,i — m c as the force acting on the interface. The 
left term of eq. (Q) is a discrete time derivative while 
+ hi-i — 2hi is the discrete Laplacian in one dimen- 
sion. After coarse-graining one obtains the continuum 
equation p6| 

Xd t h = T\7 2 h-m c + 7](x,h). (5) 

Here A and T are coarse- graining parameters, A can be 
interpreted as a viscosity coefficient and T as a surface 
tension. 

r)(x, h) is a quenched noise which may have different 
origins. For instance, the critical slope may be a random 
variable reflecting the randomness in the local rearrange- 
ments of sand grains after each toppling. In a more real- 
istic model m c is thus a random variable which changes 
its value from site to site after each toppling event. When 
h(x,t) advances, i.e. the site x topples, a new random 
slope is assigned. The randomness in the critical slope 
is thus reflected in eq. (||) through the quenched noise 
r/(x,h). The quenched noise rj(x,h) will be assumed a 
Gaussian noise with zero mean and noise correlator 

(r)(x, h)rj(x', t')) = S d (x ~ x')A(t - t'), (6) 

where A(h) is a symmetric function, i.e. A(— h) = A(h). 
Depending on the boundary conditions for A(h) one can 
distinguish two cases p^-p5|]. For random disorder A(h) 



goes to zero for large h while for periodic pinning forces 
A(h) is periodic. 

The problem will be completely defined after the initial 
and boundary conditions are specified. We are interested 
in the stationary solution of eq. (||) and therefore the ini- 
tial condition is irrelevant. Grains will be added to the 
system at constant rate c at site i — 0. Hence, if we start 
with an initial configuration where zq = the number 
of grains at site i — is obtained adding the number of 
grains received from the external field (ct) , and subtract- 
ing the number of grains transferred to the right neighbor 
(hi), i.e. 

Zi = ct - hi. (7) 

Using this expression one may compute the local slope 
m, = Zi — resulting 

rrii — m c = hi — 2ho — m c . (8) 

This evolution rule together with that in eq. (||) leads to 
the coarse grained equation 

dh 

Xdth = k— — h ct — h — m c + r](x, h), (9) 
ox 

which is the boundary condition at x — 0. However, in 
the stationary state the major contribution to the right 
hand side of this equation is given by the term ct — h and, 
therefore, this boundary condition can be approximated 
by the more simple condition 

h(0,t)=ct. (10) 

On the other hand, at x = L — 1 the boundary is open 
which implies that the site x = L will never topples, i.e. 

h(L,t) = 0. (11) 

Eq. (j|) can be generalized to a d-dimensional model. 
To be more specific we consider a d-dimensional hyper- 
cube of linear size L where grains are added at face 
xii = and the boundary X|| — L is open. The equa- 
tion of motion of h(x, t) will be given by eq. (|5|), where 
V 2 will be now the d-dimensional Laplacian, with the 
boundary conditions 

h\ Xll=0 = ct, h\ xu =L = 0- (12) 

Periodic boundary conditions are assumed in the other 
directions. Eqs. (||) and (|l^) completely defines the LID 
model dragged at one end. For c — > we recover the 
equation of motion proposed by Paczuski and Boettcher 
p2| , which describes the dynamics of critical slope mod- 
els in the limit of separation of time scales and no local 
dissipation. 
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III. STATIONARY STATE 



F = F c + const.c 1/p . 



(17) 



Eq. (||) is similar to the equation of motion of LID 
models. However, in this case the interface is dragged at 
one end instead of being driven by a constant force. The 
drag at one end makes the problem asymmetric, which 
leads to a gradient in h(x,t). It is thus easier to look 
for an expansion around this gradient. With this idea in 
mind we look for a solution of the form 



h(x,t) = h (x\\,t) + y{x,t), 



(13) 



where ho(x, t) will be determined imposing the constraint 
y\x< |=o = y\x\\=L — (symmetric boundary conditions 
for y(x,t)) and introducing a constant force term in the 
equation for y(x,t). These requirements lead to the fol- 
lowing problem for ho(x,t) 



Br 2 
° X \\ 



F — m r 



(14) 



with the boundary conditions in eq. (pL2[) . The solution 
of this problem is given by 



ho(x,t) 



F — m r r, / F — m r „\ X\\ 

c -x\ - [ ct + — L 2 ) + «'/. 



2T 
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L 



(15) 



Then, substituting eq. ( [l3|) in eq. (||), with ha(x, t) given 
by eq. (pi), and taking the limits vt > (F - m c )L 2 /2T 
and xn <C we obtain the following equation for y(x,t) 



\d t y = TV 2 y + F - Ac + v(x, ct + y). 



(16) 



This equation describe the fluctuations of the interface 
profile h(x, t) around ho(xu , t), away from the boundaries 
and for very long times. The constant force F, introduced 
in eq. ( |l4[ ) , will be determined self-consistently using the 
constraint (y(x,t)} = 0. Within this approximations the 
anisotropy introduced by the boundary conditions, such 
that the resulting equation for y(x,t) is isotropic. The 
influence of the anisotropy will be considered in the next 
subsection to compute the avalanche exponents. 

Eq. ( |l6| ) is identical to the one obtained for the fluctu- 
ations of an elastic interface driven uniformly by a con- 
stant force F [2^,Q. However, in this case F is a fixed 
parameter while the interface velocity v (c in eq. (|l6|)) is 
obtained self-consistently from the equation of motion. A 
depinning transition takes place at certain critical force 
F c determined by the disorder. For F < F c the interface 
is pinned after certain finite time while for F > F c it 
moves with finite average velocity v ~ (F — F C )P . On 
the contrary, in our model the interface velocity c is the 
fixed parameter, while F is determined self-consistently 
from the equation of motion. Since c > the system will 
always be above the depinning transition, i.e. F > F c . 
Moreover, to obtain an average interface velocity c we 
should have c ~ (F — Fc) 13 and therefore 



To reach the critical state F — F c we must then fine- 
tune c to zero. In other words, the critical state will be 
obtained in the limit of separation of time scales c — > 0, 
which is usually satisfied in numerical simulations. Ac- 
cording to eq. (15|) adjusting the constant force F we are 
just adjusting the curvature of the interface profile h(x, t) 
along the xu direction. Hence, the system self-organizes 
into a stationary state where the curvature of the inter- 
face balances the pinning forces. This conclusion, which 
was already conjectured by Paczuski and Boettcher, was 
obtaining here by solving the LID model dragged at one 
end. 



A. Scaling exponents 

The fluctuations around the average profile are char- 
acterized by the pair correlation function. At the critical 
state the pair correlation function is expected to obey the 
scaling law p3pl 



([y(x,t)-y(0,0)¥)~\x\«g(t/\xn, 



(18) 



where g{x) is an scaling function and z and £ are the 
dynamic and roughness exponents, respectively. These 
exponents are obtained through a RG analysis of eq. 
( |l6f ) |p3|-p5|. The upper critical dimension is found to 
be d c = 4 and below d c the scaling exponents are given 
by 



C 



(19) 



for a random distribution of pinning forces [|23|]24]] and 

(20) 



* = 2-i^, c=o, 



for periodic distribution of pinning Eq] . The exponents z 
and C are thus identical to those obtained for the constant 
force case. 

However, other scaling exponents result different be- 
cause of the anisotropy introduced by the boundary con- 
ditions. In the case of an interface driven by a constant 
force the average interface profile above the depinning 
transition is given by vt and, therefore, there is no pref- 
erence direction in space. On the contrary, in the case 
of an interface dragged at one end the average interface 
profile is given by ho(xu,t), which is clearly non- uniform 
along the x\\ direction as one can see from eq. (fL5|). 

To analyze the influence of the anisotropy introduced 
by the boundary conditions let us analyze the dynamics 
of active sites. Let p a (xii,x±,t) the average density of 
active sites at site (2:11, a; j_) and time t, given was a site 
active at xn = 0. Here Xi as above denotes the posi- 
tion along the preferential direction and xj_ is a d — 1 
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dimensional vector denoting the position in a plane per- 
pendicular to x\\. Since the flow of grains takes place, in 
average, along the positive x\\ direction then the average 
flux of particles through the plain x\\ — L is given by 

J(L) = J dtd d - 1 x ±Pa (L,x ± ,t). (21) 

On the other hand, at the critical state c — > the average 
density of active sites should satisfy the scaling law 

Pa( V ,,t)=t^f^,^,±y (22) 

where rj is another scaling exponent and / is an scaling 
function. Then, substituting this expression in eq. ( pi] ) 
it results that 

J(L) ~ L^^- 1 . (23) 

Now, if the system is in a stationary state for each grain 
we put at x\\ =0 one grain should go out at x = L and, 
therefore, J(L) = 1. This stationary condition will be 
satisfied only if 

(l + V )z = l, (24) 

in eq. (^. A more familiar scaling relation is obtained if 
one computes the mean avalanche size, which is given by 

(s) = J dtdx\\d d ~ 1 x^p a {x\\,x^,t). (25) 

Substituting the scaling law for p a in eq. ( p2[ ) in this 
expression it results that 

(s) ~ L {1+v)z ~ L, (26) 

which is the usual scaling dependency of the mean 
avalanche size with system size in critical slope models. 
This derivation using such simple scaling arguments is 
reported here for the first time. 

Using the RG estimates for the exponents z and ( in 
eqs. ([l9]) or (^0|) and the scaling relation (s) ~Lwe can 
compute the avalanche scaling exponents. Let s be the 
avalanche size and T its duration, which are distributed 
according to P(s) and P(T), respectively. Just at the 
critical state one expect that these distributions satisfy 
the power law behavior P(s) ~ s~ Ts and P(T) ~ T~ Tt , 
where r s and r t are the avalanche distribution expo- 
nents. However, for a system of finite size a characteristic 
avalanche size s c ~ L D and duration T c ~ L z will ap- 
pear, where D is the avalanche fractal dimension. Near 
the critical state the distributions of avalanche size and 
duration will thus satisfy the scaling laws 

P(s) ~ s- T f(s/L D ), P(T) ~ T-^g(T/L z ), (27) 

where f(x) and g{x) are some cutoff functions with the 
asymptotic behaviors f(x),g(x) ~ 1 for x <C 1 and 
f(x),g(x) < 1 for x » 1. 



The exponents t s , r t , D and z are not all independent. 
Since s ~ T z / D then the condition / dsP(s) = J dTP(T) 
implies 

(r s - 1)D = (n - l)z. (28) 

Another scaling relation can be obtained taking into ac- 
count that (s) ~ L (see eq. (|^), resulting 

(2 - t s )D = 1. (29) 

Then, from eqs. ( p8| ) and (|2^) it follows that 

r s = 2-i, r t = l + ^— i. (30) 

Finally there is a scaling relation which relates the 
avalanche dimension exponent D with the roughness 
exponent £. Below the upper critical dimension the 
avalanches are compact objects |27j and therefore s ~ 
Ahr d , where A/i is the characteristic fluctuation of the 
interface width during the avalanche and r its character- 
istic linear extent in the <i-dimcnsional substrate. Then 
since Ah ~ and s ^ r D one obtains p7| 

D = d + C (31) 

Above the upper critical dimension the avalanches are no 
more compact and D = d c = 4 [p7| . 

IV. COMPARISON WITH EXPERIMENTS AND 
NUMERICAL SIMULATIONS 

Using the values of z and C, obtained from the RG anal- 
ysis in eq. ([l9| ) or (^o|) and the scaling relations in eqs. 
d30| ) and (|31|) we can determine all the avalanche expo- 
nents. The results in one and two dimensions are shown 
in table | and [n] for random and periodic pinning, re- 
spectively. Some numerical estimates for critical slopes 
model are also shown for comparison. 

In d = 1 we count with numerical simulations of ri- 
cepile models, which are critical slope sandpile models 
with certain randomness in the toppling rule. For in- 
stance Frette |ll]] considered a ricepile model where the 
slope threshold m c is selected at random after each top- 
pling. A modified version of this model was later used 
by Christensen et al jl6). In a somehow different model 
Nunes-Amaral and Lauritsen [fl2f considered a rice pile 
model where the toppling rule is stochastic in certain 
range of slopes m c \ < m < m c i while it is determinis- 
tic above m C 2- All these ricepile models give the same 
avalanche exponents and should belong to the random 
disorder universality class. In table Q we display the more 
accurate numerical estimates reported by Nunes-Amaral 
and Lauritsen. 

In d = 2 we count with a critical slope model in- 
troduced by Bassler and Paczuski to describe the 
avalanche dynamics in superconducting vortexpiles. In 
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their model Zi is the density of vortex at site i. The evo- 
lution rules are not exactly those described above but a 
critical slope condition was used. They also considered 
quenched random pinning forces and therefore it should 
belong to the random disorder universality class. 

More recently Cruz, Mulet and Atshuler have 
made simulations of the model introduced by Bassler 
and Paczusky but considering both random and periodic 
pinning forces. In the case of random pinning their nu- 
merical estimates reproduced, within the numerical error, 
those reported by Bassler and Paczuski |15|. However, in 
the case of periodic pinning they obtained very different 
exponents (see table ||) suggesting that there are two 
universality class. Within our analysis the existence of 
these two universality classes is clear, corresponding with 
random and periodic distributions of pinning forces. 

In all cases we observe a very good agreement between 
the numerical estimates and our predictions. This agree- 
ment is more surprising because the dimensions consid- 
ered are far from the upper critical dimension d c = 4. 
In table Q we also display the scaling exponents D and 
z obtained from numerical simulations of the LID model 
driven at constant force, which are expected to be identi- 
cal to those for the LID model driven at one end. As one 
can see the agreement is quite good, event better than 
with the RG estimates. 

Recently Tadic and Nowak |^9| has observed that 
the scaling exponents of the random-field Ising model 
(RFIM) in the low disorder regime are very close to those 
of the ricepile model (see table |). In particular they 
considered a diluted two-dimensional Ising model with 
weak random fields and with an anisotropic initial con- 
dition. In the low disorder regime a single domain wall 
separating two regions with different spin orientations is 
observed. In this regime the Barkhausen avalanches are 
attributed to the fluctuations in the domain wall. The 
agreement with the exponents for the ricepile model sug- 
gests that the domain walls in this regime may be also de- 
scribed by the equation of motion ||, while the anisotropy 
was introduced through the initial conditions. On the 
contrary, on the high disorder regime a multi-domain 
structure is obtained, resulting in different avalanche ex- 
ponents. 

In this analysis we have not included a comparison with 
some available experimental results for sandpiles ri- 
cepiles [[| and vortexpiles j| because it is very difficult to 
compare our results with experimental measurements. In 
most of the experiments the measured magnitude is the 
number of grains (vortexes) leaving (entering) the system 
which will be denoted by sq. In general sq can be also de- 
scribed by the scaling law P(s ) ~ Sq T ° f(s / L D °). How- 
ever, in average, for one grain (vortex) entering the sys- 
tem one grain (vortex) goes out and therefore (sq) = 1. 
Then if one computes (sq) using the scaling law for -P(so) 
one obtains (sq) ~ £( 2 - r o)-D ^ ^ which immediately 
gives To = 2. This value is near the one reported for real 
sandpiles, ricepiles and vortex piles. On the contrary, 
in numerical simulations and in the model analyzed in 



this work the avalanche size s is given by the number of 
toppling events which is not necessarily proportional to 
the flow of grains (vortexes), the magnitude measured in 
experiments. In fact, in the case of the avalanche size s, 
the conservation law yields t s = 2 — 1/D < 2) (see eq. 
(E9|). Hence, we conclude that tq ^ t s . 



V. CONCLUSIONS 

We conclude that ricepile models, the vortexpile model 
by Bassler and Paczuski, and the RFIM in the low dis- 
order regime belong to the same universality class, that 
of critical slope sandpile models. The prototype model 
for this universality class is the LID model dragged at 
one end, as it was already conjectured by Paczuski and 
Boettcher 
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Our analytical treatment has demon- 
strated that their guess was correct. For the first time, we 
have solved the LID model dragged at one end, obtain- 
ing the complete set of scaling exponents. Some of these 
results were supported by previous RG calculations de- 
veloped for constant force LID but others, like the linear 
dependence of the susceptibility with system size, are in- 
trinsic to this model. Moreover, as in the constant force 
case, we found two different universality classes corre- 
sponding with random and periodic pinning forces. Our 
predictions were found in very good agreement with nu- 
merical simulations of different critical slope models. 
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d Model 



1 RP 
RFIM 
LID 
RG 

2 VP 
LID 
RG 



1.53(5) 
1.58 



=1.5 



1.63(2) 



n 


D 


z Ref 


1.84(5) 
1.89 


2.20(5) 
2.23 


1.40(5) 
1.45 


29 


|=1.75 


2.25 
2 


1.42(3) 
| «1.33 


26 
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TABLE I. Avalanche scaling exponents of the ricepile 
model (RP), the vortexpile model (VP), the random-field 
Ising model (RFIM), and the LID model (LID), together with 
our estimates using previous RG calculations. Error bars are 
indicated between parenthesis. 



d 


Model 


T 


n 


D 


Ref. 


1 


RG 


1 


l 


1 


1 


2 


VP 
RG 


1.45(2) 
|=1.5 


1.70(8) 
3 =1.75 


2.2(1) 
2 


1.6(1) m 

| «1.33 



TABLE II. Avalanche scaling exponents of the vortexpile 
(VP) model with periodic pinning, together with our esti- 
mates using previous RG calculations. Error bars are indi- 
cated between parenthesis. 



G 



